FETMA: A Simple Code for Ther mo-Mechanical 
Analysis on BWR Fuel Rods 

Hector Hernandez- Lopez 

Nuclear System Department, National Institute for Nuclear Research 
Carr. Mexico -Toluca s/n (km 36.5), La marques a, Ocoyoacac, Mexico, 52750 
hector.hernandez@ in in .gob . mx 



Abstract- The new designs of fuel elements are enhanced to 
reach high burn-up. If these new designs can satisfy the safety 
design constraints, a strong positive economic impact can be 
achieved. Thus, for safety and economic evaluations, the 
behaviour of fuel elements needs to be analysed. A 
computational tool was developed to analyse the thermo- 
mechanical behaviour of current new design fuel elements 
during normal or transient operation. This computational code 
solves the diffusion equation with six neutron energy groups, so 
the influence of fast neutrons on the cladding can be better 
estimated. Fuel depletion, not considering the gadolinium as 
burnable poison, calculations for fission product yield is also 
performed. A description of the code and result obtained are 
presented in this task. 
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I. INTRODUCTION 

In the nuclear industry, the new fuel assembly designs 
are continuously under development. These designs are 
submitted to high burn -up and long residence time inside the 
core of the nuclear reactor, as is the case of Laguna Verde 
Nuclear Power Plant (LVNPP). The new designs should 
provide warranty to satisfy high safety levels, given that a 
fuel failure plays an important role in safety when the 
economic impact is too high. The nuclear fuel manufacturer 
evaluates their designs to warrant the safety of the Nuclear 
Power Plant. However, due to safety and economic topics, 
the behaviour of fuel elements needs to be analyzed during 
its life inside the reactor core. In the National Institute for 
Nuclear Research (ININ, Instituto Nacional de 
Investigaciones Nucleares), a computational tool was 
developed to analyse the thermo-mechanical behaviour of 
current new designs of fuel elements during normal or 
transient operation to be used at LVNPP. 

The aim of this task is to develop a computer code for 
fuel elements thermo-mechanical analysis for assessing the 
fuel elements performance normal and transient operation. 
This code considers that different phenomena will occur at a 
fuel element when introduced into a power reactor. These 
phenomena are due to fission, a chemical and physical 
process. The processes of heat and mass transport that occur 
as a result of heat are generated by fission, and chemical 
reactions that occur between the materials are from fuel, 
cladding and coolant. 

The code is divided into three parts: neutronic, thermal 
and mechanical. From a difference respect to other codes, 
the FETMA code evaluates all aspects involved in fuel 
element performance inside reactor core, and takes into 
account several operation conditions at present. The 



sequence of calculation employed by the FETMA code is 
shown in Fig. 1. 

The neutron microscopic cross section library for six 
energy neutron groups contains 19 isotopes including 
material components of fuel rod and 4 main fission products. 
The one dimension, axial direction, diffusion equation in six 
energy groups is approximated considering fuel 
homogenization by node. The fuel rod power calculation is 
made taking into account its position in the fuel assembly 
and location in the reactor core. Also, burn-up and 
production of fission products are calculated. 




Fig. 1 Sequence of calculation employed by FETMA code 



The first approach employs neutron flux from the axial 
power distribution obtained by CM-PRESTO [1] for the 
assembly where the fuel rod is located, CM -PRESTO 
ensemble average values reported by the additional use of 
HELIOS code [2] to obtain the radial distribution in fuel 
lattice. HELIOS takes into account the presence or absence 
of the control rod and the voids ratio at which the fuel lattice 
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is considered to produce a search that implicitly takes into 
account these considerations. 

The axial temperature distribution in coolant channel 
and one or two phases along the channel are approximate. In 
this case it is necessary to execute calculations about void 
and steam quality distributions, as well as pressure drop in 
the coolant channel. The calculations for radial temperature 
distribution in the fuel rod are executed considering fuel 
density changes and gas volume is increased into the gap. 

Finally, the stresses and strains distribution is caused by 
temperature changes and materials that affect both the fuel 
and cladding, causing the phenomena of densification and 
creep. The stresses and strains distributions is caused purely 
by mechanical phenomena because of pellet-cladding 
interaction and changes in the internal pressure due to 
fission gas released from the pellet to gap. By assessment, 
the effects are caused by radiation to fuel, which is the result 
of swelling suffered by the pellets due to accumulation of 
fission gas, as well as cladding, which is presented as the 
creep of the material product of the neutron irradiation. 

II. NEUTRONIC ANALYSIS 

A. Microscopic Cross Sections 

To generate the microscopic cross section library, the 
NJOY code [3] was used and the database ENDF-B/IV-rel 4. 
[4], the library generation are obtained by the application of 
the modules contained in the NJOY code, as shown in Fig. 2. 



TABLE I ENERGY GROUP SAND LETHARGY USED IN FETMA CODE 
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Fig. 2 Sequence used in NJOY for calculation of microscopic cross sections 

The library is obtained after the reconstruction of 
parameters of resonance, characterization of the Doppler 
widening at high temperatures and energy losses, 
calculation of the self-shielding in the region of unresolved 
resonances, generation of the dispersion effective sections in 
the thermal energy interval, and finally the processing of the 
sections by energy group. 

Table 1 shows the energy groups considered in FETMA. 
These energy cuts are taken from Stammler [5]. Similarly, 
Table II shows the isotopes included in the library of cross 
sections for FETMA. As it can be seen, it includes data for 
the coolant, cladding, fission products, poisons, burnable 
poisons, fuel and plutonium 
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TABLE II MICROSCOPIC CROSS SECTION ISOTOPES LIBRARY CONTAINED IN 
FETMA CODE 
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B. Neutron Diffusion 

FETMA resolves the diffusion equation in a system of 
six differential equations coupled with six unknowns; 
therefore, to make the spatial discretization we are using, the 
cash method of integration is based on finite differences. In 
order to solve the system, it applies the Gauss-Seidel 
iterative method, as shown in the following Equation (1). 



(1) 
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C. Fuel Burn-up 

During the operation of a nuclear reactor for fuel 
material composition, it will change as the fissile isotopes 
are consumed and fission products are produced. This 
process can be monitored over the lifetime of the reactor in 
an effort to determine the material composition and 
reactivity as a function of the energy released. This requires 
the depletion and production chains for the main isotopes, 

235 2 38 

U and U, coupled with the equations that determine the 
neutron flux in the reactor core. On the other hand, the 
schemes that continuously estimate the content of fission 
products such as 135 Xe, 149 Pm and 149 Sm observed the 
importance of 135 I, 135 Xe for both neutron and for the 
mechanical performance. A good approximation is needed 
for the volume that would come to occupy within the fuel 
rod. 

Both fission products and isotopes that make nuclear 
fuel follow the pattern of generation-consumption are 
mathematically detailed in Equation (2). 



dN A 
dt 



/V 4 +AJV„ + 



N r 



(2) 



Where: Na, b, c is the isotopic density A, B and C 
isotopes. <ja,c is the microscopic cross section of absorption 
of isotopes A and C. is the decay constant of isotope A 
and B. 
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HI. TEMPERATURES DISTRIBUTION 

A. One Phase 

In boiling water reactors (BWR), the coolant 
temperature, along with the fuel rod should be studied in 
two physical regions. The first occurs when water maintains 
its liquid phase: hz<hsat. The change of temperature along 
the water channel considers the enthalpy as a function of 
system pressure and it should be based on drop pressure due 
to friction in the channel, pressure hydrostatics, and the 
form of the channel. 

The coolant's density and temperature are determined as 
a function of enthalpy (h) and pressure (p). In the region 
where coolant is boiling, (hi>h SAT or Zi>z SBC ), the 
temperature is found in the balance equation. The pressure 
drop is modified due to two phases flow effects. Changes in 
water temperature along the cooling channel are given by 
Equation (3). 



■■h, + — {q,+q M ) 
2w 



(3) 



where: q is the rate of heat generation [W]. w is the mass 
flow of refrigerant [kg / sec]. Az is the difference in height 
of the node [m]. hi is the enthalpy in the ith node [kJ / kg] . 

Where the enthalpy h is a function of system pressure at 
z, therefore variations in system pressure due to friction APf, 
the hydrostatic APh and APF channel shape must be 
considered; system pressure at z is defined by Equation (4) 
[6]. 
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where: P; is the pressure in the coolant channel at point i 
[Bar]. Dh is the equivalent diameter of the fuel rod analyzed 
[dimensionless]. f(R e ) is the fanning friction factor as a 
function of Reynolds number and relative roughness [-]. Re 
is the Reynolds number, calculated for a fluid with forced 
circulation [-]. p ; is the density of the refrigerant at point i 
[Kg/m3]. Ki is a friction factor determined for each type of 
grid spacer, g is the gravitational constant [9.8 m/s2]. u z 2 is 
the square of the average velocity of the refrigerant at point i 
[m/s]. 

Once you have determined the enthalpy, it defines the 
density and temperature of the coolant as a function of h and 
P using the correlations published by Jordan [7]. 

B. Two Phases Flow 

This is a regional characteristic of reactors cooled by 
boiling light water (BWR), which presents the boiling of 
coolant, hj = h sat where, in order to find the temperature of 
the system, is part of the energy balance equation, which 
gives us the enthalpy at each point of the channel. The 
pressure drop is modified according to the effects of two- 
phase flows; these changes are implemented into the 
Martinelli-Nelson correlation [8] for the pressure drop due 
to friction, see Equation (5). 



where: AP; f is the pressure drop in one phase, at the point z, 
Rmn is the Martinelli-Nelson multiplier, depending on the 
pressure and the quality of the coolant at point Zj. 

Moreover, the pressure drop by the acceleration APa has 
to be taken account of and has the following form, see 
Equation (6). 



V * p,u 



Az 



(6) 



where: a is the fraction of voids in the coolant, pi is the 
density of the refrigerant in its liquid phase. p g is the density 
of the refrigerant in its gaseous phase 

Once the system pressure is found, the correlations are 
applied; it is developed by Saul [9] to find h g , hi, pi, p g . 

In addition to T sat Tj = (Pi). 
C. Radial Temperature Distribution 

The energy transport equation describes the temperature 
distribution in a solid, which is assumed to be 
incompressible and with negligible thermal expansion. Also 
a steady state is assumed so that the resolved expression for 
the temperature distribution in the radial direction is given 
by Equation (7). 



-Vk(r,T)VT(r)=q'"(f) 



(7) 



(5) 



where: k(r,r) is the coefficient of thermal conductivity as a 
function of position and temperature. 

For purposes calculation, one can consider the surface 
temperature of the cladding as a reference temperature; it 
will be determined in terms of coolant temperature. To 
simplify the calculations, certain assumptions were made. 

We assume that one can neglect thermal conduction in 
the axial direction; this is because temperature gradients in 
the radial direction are greater by several orders of 
magnitude to gradients in the axial direction. As the heat 
transfer in the axial direction will take place outside of the 
fuel by the coolant through forced circulation. The fission 
energy appears as a uniform heat source distributed through 
the combustible material; usually there is some variation in 
the heat source because it is proportional to the neutron flux 
in the fuel. Only the heat transfer is considered in steady 
state. 

Given these assumptions, for the implementation of 
Equation (7) in the temperature analysis, the different 
coefficients of thermal conductivity must consider, as 
presented in the fuel materials, i.e, the thermal conductivity 
of the pellet fuel, thermal conductivity, which occurs in the 
pellet -cladding gap and the thermal conductivity through the 
cladding. 

FETMA code contemplates the use of four models to 
represent the coefficient of thermal conductivity, which are 
functions of density, burn -up and amount of gadolinium. 
The first model to calculate the thermal conductivity is 
taken from the properties of nuclear materials covered 
MATPRO Version 9 [10]. The second model is determined 
from experimental data obtained in the reactor HALDEN 
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[11]. The third model includes dependence on the density 
and temperature of the fuel, which was developed by Lucuta 
[12]. The final model considers a dependency on the content 
of gadolinium in the fuel pellet developed by Fukushima 
[13]. Fig. 3 shows the behavior of models of U02 thermal 
conductivity used in FETMA code. 




Temperature |"C) 



Fig. 3 Models of U02 thermal conductivity used in FETMA code 

D. Thermal Conductance Gap 

Since some fuel rods use He or other gas to fill the 
clearance between combustible material and the cladding, 
the heat transfer model through it should consider the 
conductance through the gas and thermal radiation between 
surfaces, combustible material and cladding. In addition to 
heating and irradiation, swelling occurs in fuel material, 
which can produce contacts between fuel and cladding. 

Since the problem of determining the conductance in the 
clearance (hgap) is very difficult to resolve by the number 
of variables that are involved in design, hgap is considered 1 
W/m2°K. However, there are relationships which may give 
a more realistic approximation of the value of hgap; in our 
case we take a model proposed by Ross and Stoute [14]. 

E. Cladding 

To calculate the temperature of the cladding, the 
coefficient of thermal conductivity of zircaloy is considered 
so that Anderson [15] evaluates test data and adjust 
minimum through a cubic equation in terms of temperature. 
Equation (8) can be applied both to zircaloy-2 and zircaloy- 
4, as the fit of the data was carried out with data from both 
alloys. 

k(T) = 7. 151 + 2.472x1 OT + 1. 674x1 OT 2 - 3.334x1 0""T 3 (8) 

where: T is the temperature in the cladding [°K]. K is the 
coefficient of thermal conductivity [W/m°K] . 

IV. THERMAL-MECHANICAL STRESS 
A. Fuel Strain 

In the case of the ceramic fuel, the problems of thermal 
stress are of primary importance. In general, ceramics have 
little ductibility compared with metallic materials, and 
therefore are more fragile. Thus, the fuel rods are subject to 
another type of stress, which is produced by swelling in the 
pellet and material creep. This stress should be considered 
in the mechanical analysis of the fuel rod. 



Equation (9) shows the model used for calculating the 
deformation e in terms of stress a, in cylindrical coordinates. 

= ^ k - u k + ct z )] + + 6„ + ** 
e e =— k - y k + °" z )] + s „ + c + £ ,„ (9) 
e, = £ k - u k + °v )] + k + £ .„ + 

where: £ is the strain in the direction i (r, 9, z) [m / m]. a; 
stress in the direction i (r, 9, z) [MPa]. E is Young's modulus 
[MPa]. v is the Poison's modulus [dimensionless]. e cr is 

caused by creep deformation [ml m]. s m \& the strain caused 

by swelling [m / m]. e iS is the deformation due to thermal 

expansion, given by aT [ml m]. A is the thermal expansion 
coefficient [ml (m° C)]. 

It should be noted that strain due to creep, swelling and 
thermal expansion does not depend on the position and can 
be considered constant. To solve this we must assume that 
the axial deformation of the fuel is absorbed by the springs 
that hold the fuel inside the bar. 

Cladding of the thermal stress is due to an expansion 
produced by a differential temperature between the inner 
and external walls of the rod. The region inside the cladding, 
near the fuel, is hotter than the external region, which causes 
compression stress near the inside face and tensile stress 
near the external face of the cladding. 

The stresses, due to pressure, are the result of the 
pressure difference between inside and outside of the fuel 
rod. The external pressure is caused by the pressure of 
system operation, while the internal pressure is caused by 
the expansion of the gas -filled rod and the release of fission 
gases as a result of irradiation. 

B. Sweeling 

During the same process of burning a fuel, each fission 
atom is replaced by two fission product atoms, with about 
half the atomic weight of the fissionable atom each. Of these 
atoms, about 15% are inert gases, such as Xe and Kr, and 
the remaining are solids and soluble gaseous oxygen. 

Once the initial densification is completed, the fuel tends 
to swell as the burn-up is increasing. This phase of swelling 
is sometimes less pronounced at lower densities, i.e, a 
greater number of pores in the material versus a combustible 
material with fewer pores, which represents a higher density. 

At specific temperatures and operation levels due to the 
degree of burn-up of a light water reactor, the fuel-swelling 
ratio has a value of approximately 0.5% of volume per 10 20 
fissions per cubic centimeter of fuel material, approximately 
4.15 GWD/TMU (10.GWD/TMU= 2.41xl0 20 fissions/cm 3 ). 

The volume change caused by non-gaseous fission 
products is very difficult to measure; however, there had 
been many experiments to establish relative amounts, taking 
into account the elements and compounds produced. The 
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code makes use of the model described by Equation (10), 
proposed by [16, 17]. 



= 7.435x1 O^B.-B,,) 



(10) 



where: Bj is burn-up in a time step i [MW/KgU]. P is the 



initial density of fuel material [kg/m 



AV~\ 



is the fraction 



of change in volume due to solid fission products. 

The swelling due to fission gases is caused by the 
increased number of bubbles of gas produced; however, the 
physical mechanisms that cause it are too complex. As in 
the swelling caused by solids, gaseous products are 
calculated by an experimental correlation; Equation (11) 
expressed in terms of temperature and burn-up of the fuel 
pellet [18, 19]. 



= 2.617x10-" pB(2800 - r)exp[- 0.0162(2800 - r)]exp(- 2.4x1 0" pB) (ID 



where: T is the average fuel temperature [°K]. B is burned 
in a time step [MW/KgU]. P is the initial density of fuel 

material [kg/m3]. ( is the fraction of change in volume 

due to gaseous fission products. 

Therefore, the swelling caused by the pellet fuel fission 
products is the sum of both cases [20]. 

On the other hand, the swelling that occurred in the 
cladding is approximated by the model used by Benjamin 
[21] in the code ISUNE2 and is described by Equation (12). 



AV 



= 9x10 "(0-t)' 



(12) 



(4.028-3.712x10 2 7/ + 1.0145xl0 'T 1 -7.879x10 ! 7" 



where: T is the average temperature in cladding [°C]. T is 
the irradiation time [seg]. </> is the fast neutron flux with 
energies > 1 MeV. [n/cm 2 -seg] 



C. Fission Gas Release 



As mentioned above, as a result of the 235 U fission, the 
Xe is produced, an insoluble inert gas in the fuel matrix can 
be contained within the fuel matrix causing swelling of the 
matrix, or else be released to the cavities forming the fuel 
cladding gap and the central void of the tablet, if any, 
through the cracks that are formed within the fuel element 
and connecting pores too large, with the above mentioned 
regions. 

However, only a fraction of these Xe atoms produced, 
manage to escape to these areas, thereby producing an 
increase in filling gas pressure that is within the recoil, 
which has a direct impact on the tensions experienced by the 
fuel cladding. 

The fractions of gas release are very difficult to measure; 
however, several experimental measurements support that 
these quantities are evaluated by the model proposed by [22]. 

D. Creep 

Creep is the tendency of a solid material to slowly move 
or deform permanently under the influence of stresses. It 



occurs as a result of long term exposure to levels of stress 
that are below the yield strength of the material. Creep is 
more severe in materials that are subjected to heat for long 
periods and near the melting point. Creep always increases 
with temperature. 

The rate of this deformation is a function of the material 
properties, exposure time, exposure temperature and the 
applied structural load. Depending on the magnitude of the 
applied stress and its duration, the deformation may become 
so large that a component can no longer perform its function. 
Creep is usually of concern to engineers and metallurgists 
when evaluating components that operate under high 
stresses or high temperatures. Creep is a deformation 
mechanism that may or may not constitute a failure mode. 
Moderate creep in concrete is sometimes welcomed because 
it relieves tensile stresses that might otherwise lead to 
cracking. 

Unlike brittle fracture, creep deformation does not occur 
suddenly upon the application of stress. Instead, strain 
accumulates as a result of long-term stress. Creep 
deformation is a "time -dependent" deformation. 

In some materials, especially polycrystalline like U0 2 , 
creep approaches a continuous distribution of modes of 
sliding points of the lattice dislocation. In this case, 
incorporating the model more accurately represents the 
phenomenon of creeping, which was originally proposed by 
Bohaboy [23] and modified by Solomon [24] as shown in 
Equation (13), and it is actually used for other codes . 



(A, + A 2 F)a,_ expj^- ^ j A t aZ exp|^- - R j 
(A, +£>)G 2 + A„ + D 



(13) 



exp - 



ft 



RT 

Where: 
A,= 3.919xl0 5 
A 2 =1.305xl0" 13 
A,= -87.7 
A 4 = 203.7 
A 6 =-90.5 
A 7 = 3.723x10 



29 



Qj= 376,740 (J/mol) 

g 2 =552,552 (J/mol) 

Qj= 21,767 (J/mol) 

eWs-165.5 O 

R= 8.314 (J/mol°K) 

Tis the temperature (°K) 

D is the density fraction between [0.92, 0.98] 

G is the grain size ~ 30 fim 

The cladding creep model that was proposed by Ibrahim 
[25] is given by Equation (14). 



5.7xlO- 12 (l + 1.0xlO s r- 7 «J°" 5 ^o- tf/ +7 



lO.Oe 



(14) 
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where: T is the irradiation time [seg]. ^4 is the fast neutron 
flux with energies > 1 MeV [l/cm 2 -seg]. T is the average 
temperature in the cladding [°K]. cr, (( is the effective stress in 
the cladding [Pa]. R is the ideal gases constant. 

V. VALIDATION OF CODE 

In order to verify the results provided by code FETMA, 
an exercise was made by which the corresponding data to 
the distribution of the neutron flux were observed, as well as 
the fuel center temperatures, the fraction of fission gases 
released to gap and the elongation of the cladding. In Table 
HI, the data used for the verification exercise can be found. 

TABLE IIIPARA METERS USED IN EXERCISE OF FETMA CODE VERIFICATION 



Enrichment 


3.38 7o 


Pellet Diameter 


1 .0439 cm 


Internal Cladding Diameter 


1.0650 cm 


External CladdingDiameter 


1 .2250 cm 


Pitch 


1 .6256 cm 


Density 


96.5% 


Size grain 


7.83 urn 


Feedwater Temperature 


286 °C 


System Pressure 


70 Bar 


Coolant Flow 


1400kg/(m 2 s) 



Fig. 4 shows the neutron flux behavior for energy group 
higher than 0.821 MeV, whereas Fig. 5 shows the neutron 
flux behavior for energy group less than 0.625 eV, being 
observed that the flux represents the behavior of the power 
distribution. 




Fig. 4 Neutron flux for energy more than 0.821 MeV evaluated for 
FETMA 




To verify the capability of thermal analysis of FETMA, 
code calculation has been compared with the data of fuel 
center temperature within experiment IFA-597.2 [26]. The 
data set covers the temperature vs. power for Rod 8 during 
the first 4 ramps. Fig. 6 shows a comparison between the 
FETMA calculations and measurements of fuel center 
temperature. Fig. 7 shows the fuel center temperature 
behavior as function of the linear heat rate. 
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Fig. 6 Comparison between measured and calculated temperatures of fuel 
centerline 
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Fig. 7 Temperature at fuel centerline during power increase 

To verify the capability of the FETMA code and to 
calculate the fission gas release, the fission gas release 
behavior during irradiation has been observed in Fig. 8 and 
it shows good agreement with the theory. 

In order to verify the capability of global mechanical 
analysis of FETMA, the calculations of cladding elongation 
during first four ramps is shown in Fig. 9. 
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Fig. 5 Neutron flux of energy less than 0.625 eV evaluated for FETMA 
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Fig. 9 Cladding strain as function from maximum power duringpower 
increase 

VI. CONCLUSIONS 

A fuel behavior model in code FEIMA has been 
developed for the analysis of neutronic, thermal and 
mechanical behavior of boiling water fuel rod during 
steady-state and transient conditions. The capability of 
FETMA has been verified using some data base obtained in 
the NEA Data Bank. The results of comparison between 
FETMA calculations and experimental data for fuel center 
temperature showed that the agreement was satisfactory. 
Other parameters were just evaluated theoretically. However, 
the results obtained with FETMA were consistent with 
theoretical results. 
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